Pharmacodynamic Modeling of Acute and Chronic Effects of Methylprednisolone on Hepatic Urea Cycle Genes in Rats

Corticosteroids (CS) regulate many enzymes at both mRNA and protein levels. This study used microarrays to broadly assess regulation of various genes related to the greater urea cycle and employs pharmacokinetic/pharmacodynamic (PK/PD) modeling to quantitatively analyze and compare the temporal profiles of these genes during acute and chronic exposure to methylprednisolone (MPL). One group of adrenalectomized male Wistar rats received an intravenous bolus dose (50 mg/kg) of MPL, whereas a second group received MPL by a subcutaneous infusion (Alzet osmotic pumps) at a rate of 0.3 mg/kg/hr for seven days. The rats were sacrificed at various time points over 72 hours (acute) or 168 hours (chronic) and livers were harvested. Total RNA was extracted and Affymetrix® gene chips (RG_U34A for acute and RAE 230A for chronic) were used to identify genes regulated by CS. Besides five primary urea cycle enzymes, many other genes related to the urea cycle showed substantial changes in mRNA expression. Some genes that were simply up- or down-regulated after acute MPL showed complex biphasic patterns upon chronic infusion indicating involvement of secondary regulation. For the simplest patterns, indirect response models were used to describe the nuclear steroid-bound receptor mediated increase or decrease in gene transcription (e.g. tyrosine aminotransferase, glucocorticoid receptor). For the biphasic profiles, involvement of a secondary biosignal was assumed (e.g. ornithine decarboxylase, CCAAT/enhancer binding protein) and more complex models were derived. Microarrays were used successfully to explore CS effects on various urea cycle enzyme genes. PD models presented in this report describe testable hypotheses regarding molecular mechanisms and quantitatively characterize the direct or indirect regulation of various genes by CS.


Introduction
Glucocorticoids (GC) are a class of steroid hormones produced by the adrenal glands that have diverse effects on carbohydrate, lipid, and protein metabolism. While they are important regulators of many developmental processes, immune system and tissue integrity (i.e. bone homeostasis), one of their major physiological roles is the regulation of systemic fuel supplies. This includes gluconeogenesis, the production of glucose from non-carbohydrate sources (primarily glycerol from fat and amino acids from protein). When blood glucose from dietary sources and from breakdown of liver glycogen becomes inadequate, it is the gluconeogenic process that maintains glucose homeostasis. In order to meet this need, GC output occurs in a circadian pattern that mirrors the daily species requirement for increased glucose production. In situations of stress (where energy demands are increased), the adrenal output of GC is increased above normal levels in order to meet the demand.
The effects of GC on the immune system are exploited clinically in the use of synthetic corticosteroids (CS) as immunosuppressive/anti-infl ammatory agents. Corticosteroids are widely used to treat a variety of conditions, including organ transplantation, rheumatoid arthritis, lupus, and irritable bowel syndrome. However, such long-term use is associated with serious adverse effects that often refl ect a magnifi cation of physiological hormone actions. The development of hypertension, steroid-induced diabetes, atherosclerosis, and muscle weakness can result from chronic CS therapy and to a large extent limits their usefulness as clinical agents. Ascertaining the complex multi-gene, multi-tissue effects of steroid treatment will provide a better understanding of how pharmacologically magnified steroid-induced alterations in gene expression result in complex systemic pathologies.
Corticosteroid-induced gluconeogenesis and insulin resistance provide the foundation for steroid diabetes. Only two tissues, the liver and kidney, are capable of gluconeogenesis. Gluconeogenesis from amino acid carbon in the kidney requires disposal of the resultant ammonia as ammonium ion. However, in the liver, gluconeogenesis from amino acid carbon requires detoxification of ammonia. This is accomplished in the liver by the urea cycle, which converts ammonia to urea. The general schematic of the urea cycle, adopted from Voet et al. (Voet and Voet, 1990), is shown in Figure 1. Urea synthesis alone requires five enzymes (arginosuccinate lyase, arginosuccinate synthetase, arginase, ornithine transcarbamoylase, and carbamoyl phosphate synthetase). In addition to these fi ve genes, we have defi ned the "greater urea cycle" and also examined the expression of mRNA for nine additional enzymes that play a role in the disposal of nitrogen in the liver (glutamine synthetase, glutamine dehydrogenase, malate dehydrogenase, carbonic anhydrase, glutaminase, ornithine decarboxylase, tyrosine aminotransferase, aspartate aminotransferase and alanine aminotransferase), and two additional genes involved in transcriptional regulation of those genes; glucocorticoid receptor (GR) and CCAAT/enhancer binding protein (CEBP-β).
Most GC effects are mediated at the molecular level by changes in the expression of specifi c mRNAs. Gene arrays can provide a method of high throughput data collection that is necessary for constructing comprehensive information on the transcriptional basis of complex systemic polygenic phenomena. When microarrays are used in a rich in vivo time series experiment they yield temporal patterns of changes in gene expression that illustrate the cascade of molecular events in time that mediate broad multi-gene responses. Previously we described the mining and analysis of microarray time series illustrating the responses of liver, skeletal muscle, and kidney taken from the same set of animals to a single bolus dose of methylprednisolone (MPL) (Almon et al. 2005a;Almon et al. 2004;Almon et al. 2005c). These time series included individual chips from multiple control animals as well as multiple animals at each of sixteen times over a 72 hour period following bolus dosing with MPL.
Because these experiments were initiated using adrenalectomized animals, the drug in essence acts as a stimulus that perturbs the homeostatic balance of the system, and the deviation of the system and its return to the original state were monitored. Among the genes identifi ed as GC regulated were the 16 genes previously designated as the "greater urea cycle". Although very useful, a single time series only views the dynamics of the system in response to the one stimulus. A pharmacological time series differs from most time series studies (for example those assessing developmental changes) in that it can be evaluated using a different dosing regimen. The results from both phases can be used to group genes into clusters with common mechanisms of regulation. If two or more genes have a common mechanism of regulation, then their response profi les should be the same regardless of the dosing regimen. With such data, the challenge then becomes one of constructing rational, quantitative, mechanism-based frameworks that describe the relationships between the elements of the cascades. Kinetic/dynamic modeling provides quantitative, testable, mechanism-based hypotheses concerning the relationship between drug kinetics and elements of cascades that continue long after the drug has dissipated. Such models can accommodate hierarchal cascades where one process generates effectors or mediators for other processes. They can also accommodate convergence of cascades that commonly occur in the control of the expression of genes where binding sites for multiple nuclear factors participate in the regulation of the level of expression of a particular mRNA. This mRNA in turn becomes an endogenous mediator for the expression of proteins that may become effector molecules for other processes. Here, we present the dynamic picture of steroid modulation of expression of the "greater urea cycle" (Fig. 1).

Experimental
Livers were obtained from animal studies performed in our laboratory (Ramakrishnan et al. 2002;Sun et al. 1998a). In the acute study, male ADX Wistar rats weighing 225-250 g (Harlan Sprague-Dawley Inc., Indianapolis) were subjected to right external jugular vein cannulation under light ether anesthesia one day prior to the study. A single intravenous (IV) bolus dose of 50 mg/kg MPL (Solu-Medrol, Pharmacia-Upjohn Company, MI) was given to 47 animals via the cannula. Rats were sacrifi ced by aortic exsanguination at 0. 25, 0.5, 0.75, 1, 2, 4, 5, 5.5, 6, 7, 8, 12, 18, 30, 48 and 72 hr. The sampling time points were selected based on previous studies of receptor dynamics and enzyme induction in muscle and liver (Ramakrishnan et al. 2002;Sun et al. 1998a;Sun et al. 1998b). Four vehicle-treated rats were designated as controls and were considered as sacrifi ced at time point zero.
In the chronic-dosing study, 40 male ADX Wistar rats weighing 300-350 g were given 0.3 mg/kg/hr infusions of MPL reconstituted in supplied diluent. The infusions were given using Alzet osmotic pumps (Model 2001, fl ow rate 1 ul/hr; Alza, Palo Alto, CA). The pump drug solutions were prepared for each rat based on its pre-dose body weight. Pumps were implanted subcutaneously between the shoulder blades on the back. Four rats were sacrifi ced by aortic exsanguination at 6, 10, 13, 18, 24, 36, 48, 72, 96, and 168 hr. Four control rats were implanted with a saline-filled pump and sacrificed at selected times throughout the 7-day study period.
Blood samples were collected from the abdominal aorta, centrifuged to harvest plasma, and stored at −80 °C until analysis of plasma MPL. Livers were rapidly excised, quickly frozen in liquid nitrogen, and stored at −80 °C.
Drug assay MPL plasma concentrations were measured by a previously described normal-phase high performance liquid chromatography (HPLC) method (Haughey and Jusko, 1988). The lower limit of quantifi cation was 10 ng/ml with inter-and intra-day assay variability less than 10%.

RNA extraction and labeling
Frozen livers were ground into powder using a mortar and pestle chilled with liquid nitrogen. Total RNA extractions were carried out by a Trizolchloroform based extraction method. About 100 mg of ground liver from each animal was added to prechilled Trizol Reagent (Invitrogen, Carlsbad, CA) at a ratio of tissue: Trizol of 1:10. Extractions were performed according to manufacturer protocols. Extracted RNA was further purifi ed by passage through RNAeasy minicolumns (Qiagen, Valencia, CA). Final extracted RNA samples were resuspended in nuclease-free water. Quantity of total RNA was determined by spectrophotometry and purity was assessed by agarose gel electrophoresis. Extracted total RNA preparations were stored at −80 °C.

Microarrays
Isolated liver RNA from each individual animal was used to prepare biotinylated cRNA target according to manufacturer protocols. Then biotinylated cRNAs from the acute study were hybridized to 51 (47 treated and 4 control animals) individual Affymetrix GeneChips ® Rat Genome U34A (Affymetrix, Santa Clara, CA), which contained 8799 probe sets. The target cRNAs from the chronic study were hybridized to 44 (40 treated and 4 control animals) individual Affymetrix GeneChips ® Rat Genome 230A (Affymetrix, Santa Clara, CA), which contained 15967 probe sets. Oligonucleotide microarrays were utilized because of their high reproducibility between separate arrays. The entire data set was submitted to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus database (GDS253 and GDS972) and is also available online at http://pepr.cnmcresearch.org/ . Our approach to identifying probes of interest has been described in previously published articles on mining the datasets for muscle, liver, and kidney from the acute MPL treated animals (Almon et al. 2005a;Almon et al. 2004;Almon et al. 2005c). Literature searches of the 1512 probe sets identifi ed as regulated in the acute liver dataset yielded 16 probe sets for 16 genes relevant to ammonia detoxification. Probe sets for these 16 genes were also present on the 230A chip and used to analyze the livers following chronic infusion. These probe sets were statistically examined by ANOVA with a Tukey posthoc test (p Ͻ 0.05) using GeneSpring 7.0 software (Silicon Genetics, Redwood City, CA) and were signifi cantly different when compared with the control animals.
Data from both acute and chronic studies were simultaneously modeled. The data for each probe set were transformed so that the values for all probe sets were within the same range. To accomplish this, individual probe set values on each chip were divided by the mean of the four control values for that gene. Thus, the "normalized intensity" for each probe set has a value of 1 at zero time and either decreases, increases, or remains the same when compared to the controls over the time series.

Pharmacokinetics
The PK of MPL for both the bolus and infusion regimens were described by a two-compartment model with a zero-order input k 0 into the central plasma compartment for the infusion regimen as given below: where A p and A t are the amounts of drug in the plasma and tissue compartments. The k 12 , k 21 and V p are the distribution rate constants and central volume of distribution; CL is clearance. These parameters were fi xed based on previous literature values (Ramakrishnan et al. 2002;Sun et al. 1998a). For the i.v. bolus kinetics, the same equations were used with deletion of k 0 and use of A p (0) = Dose. The kinetic parameters were fi xed (Table 1) and were used as the driving force for the dynamic modeling.

Receptor dynamics
Glucocorticoid receptor dynamics in liver following single IV 50 mg/kg MPL dosing was previously described by a mechanistic receptor-gene mediated pharmacodynamic model (Fig. 2). The present model was developed using the assumption that being moderately lipophilic, unbound MPL diffuses through cell membranes and binds with cytosolic free glucocorticoid receptors (GR). Drug-receptor complex (DR) then translocates into the nucleus (nuclear DR: DR(N)) where it dimerizes and binds to a specifi c glucocorticoid response element (GRE) in the target DNA. The binding of DR(N) and GRE enhances or inhibits the expression of target genes. The CS are known to inhibit the expression of their own receptors (Dong et al. 1988) by homologous down-regulation. After dissociation from DNA, GR receptors are recycled into cytosol, where receptors are either degraded or further activated by MPL (Oakley and Cidlowski, 1993;Sun et al. 1998a).
The equations describing the model are: m mRNA GR (3) where symbols represent cytosolic free plasma MPL (C MPL,f ) where the unbound fraction of MPL was fi xed to 0.23 (Kong and Jusko, 1991), free glucocorticoid receptor density (R), GR mRNA (mRNA GR ), cytosolic GR-drug complex (DR), GRdrug complex in nucleus (DR(N)), zero-order GR mRNA synthesis rate constant (k syn,GRmRNA ), and fi rst-order rate constants for GR mRNA degradation (k d gr,GRmRNA ), GR synthesis (k syn,GR ) and GR degradation (k d gr,GR ). Other rate constants include: second-order association rate constant of GR and drug (k on ), translocation rate constant of DR from cytosol to nucleus (k T ), and the overall turnover rate of DR(N) (k re ). The fraction of GR that is recycled from nucleus to cytosol (R f ) can be further activated by association of drug. Thus (1−R f ) is the fraction of GR that is degraded during one cycle of drug-receptor complex translocation. The IC 50 is the concentration of DR(N) at which the synthesis rate of GR mRNA is reduced to 50% of its baseline level. Assuming that in the absence of drug the mRNA GR remains constant, the synthesis rate is: Also in the absence of drug assuming that R 0 is constant, the translational rate constant k syn,GR is: The values for the baselines of both the GR mRNA and the receptors were fi xed and were obtained from Sun et al. (Sun et al. 1998a;Sun et al. 1998b).
The baseline values for DR and DR(N) are zero. The detailed analysis of receptor dynamics (parameters given in Table 1) were done elsewhere (Hazra et al. 2007) and were fi xed to drive further dynamic modeling.

Pharmacogenomic modeling
In the present report we measured drug kinetics, GR concentration, GR mRNA concentration, and the normalized concentrations of 16 mRNAs, some of which are transcription factors that in turn become factors controlling the expression of genes.
Since mRNAs have natural turnover processes, we have used various inhibitory or stimulatory models (Dayneka et al. 1993) to describe the primary or secondary regulation of these genes by CS.
Most of the metabolic effects of CS have been attributed to its classical model of action where binding of nuclear hormone-bound GR complexes to the regulatory regions of various target genes causes stimulation or inhibition of expression of these genes. Many of the genes related to nitrogen metabolism are known to have GREs in their promoter regions and therefore presumably are directly affected by binding of the nuclear drugbound receptor to these regions. However, it has also been found that CS causes changes in gene as well as protein expression of transcription factors such as C/EBP-β or HNF-3 (Gotoh et al. 1997;Morris, 1992;Schoneveld et al. 2004) which are also strong regulators of various urea cycle genes. As a result both primary (mediated by GR binding GREs in promoter regions of target genes) as well as the secondary transcriptional regulation are passive mechanisms for GC modulation of mRNA expression.
Both mechanisms were taken into consideration to develop models to describe the dynamics of various genes related to the greater urea cycle after acute and chronic MPL. Many of the genes showing simple up-or down-regulation after acute dosing seemed to have complex patterns after chronic dosing. Our aim was to build models which were able to simultaneously describe the time profiles of these genes after both treatments.
In the absence of drug, all the genes were assumed to be synthesized at a zero-order rate (k syn,Enzm ) and degraded with a fi rst-order rate constant (k dgr,Enzm ) as follows: Since the rats were adrenalectomized, the genes under consideration were assumed to have a stationary baseline (mRNA 0 ) in the absence of the drug (control animals), thereby allowing k syn,Enzm to be calculated from: All mRNA baselines were fi xed to 1.0 except for the cases where estimation of this parameter improved the model fi ttings signifi cantly (assessed by visual inspection and AIC, SC criteria).

Stimulation of transcription
This characteristic was assumed to be mediated by the stimulatory effect of DR(N) on the synthesis rate of mRNA. Two conditions may arise; either the stimulation by DR(N) may be linear, nonsaturable ( Fig. 3A) in nature for a given gene or dose of the drug (Eq. 11a) or it could be a saturable function (Fig. 3B) where the effect of the drug is limited by capacity-limited factors such as S max and SC 50 (Eq. 11b) as given by: where S is a linear stimulation constant by which DR(N) increases the synthesis of the enzyme mRNA, SC 50 is the concentration of DR(N) responsible for 50% of maximum stimulation, and S max is the maximum stimulation capacity.

Biphasic regulation
Two different models were developed to describe biphasic regulation. In the fi rst model as depicted  in Figure 3C, DR(N) linearly stimulates the production of the enzyme mRNA whereas a biosignal (BS) derived from DR(N) is responsible for secondary inhibition of the production as given by: where k e is a fi rst-order transduction rate constant accounting for the delayed formation of BS from DR(N), S is a linear stimulation constant by which DR(N) stimulates the production rate constant, IC 50 is the concentration of BS responsible for 50% inhibition of the mRNA synthesis rate, γ is a factor describing amplifi ed effects of BS on the production of mRNA. In the absence of the drug, the BS value is set to zero. In the second model, described in Figures 3D  and 3E, the DR(N) increases the degradation of mRNA, whereas the biosignal (Eq. 12) derived from the DR(N) is responsible for a delayed increased (either linear or non-linear) production of the mRNA as given by: where S BS is the linear constant, S max and SC 50 are the capacity-limited stimulatory constants by which BS increases the production rate constant, and S DR(N) is the linear stimulation constant by which DR(N) regulates the degradation rate constant for mRNA.

Down-regulation
As depicted in Figure 3F, down-regulation of mRNA was described by: where IC 50 is the DR(N) concentration at which mRNA synthesis rate reduces to 50% of its baseline value.

Scaling Factor
Since two types of chips were used for the acute and chronic studies, a factor (SF) was incorporated in the chronic study data (Y chronic ) to account for the possible difference in sensitivity between the probe sets (Yao et al. 2007): where mRNA chronic is the data (from 230A chips) for the gene of interest after chronic infusion of MPL. Data for each gene from both the acute and chronic studies were fi tted simultaneously with ADAPT II software (D'Argenio, 1997) using the Maximum Likelihood method, where the variance model was described as: where Y represents the predicted value; σ 1 and σ 2 are the variance parameters which were fi tted, and θ represents the structural parameters. The goodness of the fi t was assessed by model convergence, visual inspection, Akaike Information Criterion (AIC), Schwarz Criterion (SC), estimator criterion value, and examination of residuals.

Pharmacokinetics and receptor dynamics
The experiments were initiated by MPL dosing. Although the kinetic data for both regimens have been previously published (Ramakrishnan et al. 2002;Sun et al. 1998a), simulated pharmacokinetic profi les for both are provided as they represent the initial driving force for the dynamic models (Fig. 4). Figure 4a shows the simulations of MPL kinetics, DR, and DR(N) following the acute dose of MPL. Although the experiment continued for 72 hours, the fi gure is truncated at 20 hours because MPL reached its lower limit of quantifi cation by 6 hr and the drug-receptor complexes returned to their baseline prior to this time. Figure 4b shows the simulation of GR and GR mRNA which is also based on previously published data (Sun et al. 1998b). Figure 4c shows the simulations of MPL kinetics, DR, and DR(N) during chronic infusion of MPL. This graph only shows the fi rst 48 hours of the 168 hour infusion period because all three variables reached a steady-state by this time. Figure 4d shows the simulation of GR and GR mRNA during chronic infusion which is also based on previously published data (Ramakrishnan et al. 2002). The steady-state concentration after the chronic dosing (i.v. infusion) was 100-fold less than the bolus dose study. The clearance (3.48 to 5.61 L/hr) difference could be attributed to the saturation of the metabolizing enzymes at the high dose. The DR(N) is used as the driving force regulating many of the genomic effects of CS.

Pharmacogenomic modeling
The dynamics of most of the genes showed more complex regulation during chronic infusion of MPL compared to acute bolus dosing. Based on the time profi les from the acute dosing study the genes of the greater urea cycle could be divided into three main classes, up-regulated, biphasic, and down-regulated. These could be further divided into various sub-categories based on the chronic dataset.

Class 1
This class is comprised of four genes which were up-regulated in response to both acute and chronic administration of MPL. Model A or B (shown in Fig. 3), where the transcription is assumed to be increased by DR(N), was used to capture the general trend of the data (Fig. 5). The dynamics of A) B)

Figure 4. Time profi les of various components in MPL pharmacokinetics and receptor dynamics.
Lines are simulations with the model given in Figure 2 after 50 mg/kg MPL i.v. injection (4A and B) or 0.3 mg/kg/hr chronic infusion for seven days (4C and D) in ADX rats using parameters listed in Table 1.
tyrosine aminotransferase (TAT) was described by Model A as discussed in previous reports based on measurements by other methods (Ramakrishnan et al. 2002;Sun et al. 1998a). The degradation rate constant (k dgr,Enz ) for TAT mRNA, 0.503 hr −1 closely resembled the reported value of 0.383-0.533 hr −1 (Ramakrishnan et al. 2002;Sun et al. 1998a). The other three genes in this class, aspartate aminotransferase (aspAT), argininosuccinate synthetase (ASS), and malate dehydrogenase (MDH), were described more accurately by Model B with capacity-limited stimulation by DR (N). The apparent tolerance observed for TAT and aspAT (Fig. 5a) was not observed for ASS or MDH (Fig. 5b), which can be explained by much lower values for SC 50 (Table 2) for the latter two genes (4.1 and 14.0 nM) compared to aspAT (107.4 nM).
The variation in the degree of stimulation (normalized signal ranging from 2.5-12) of the messages was also refl ected by their respective S max values, ranging from 0.76 to 7.3 (fmol/mg protein) −1 . The precisions of the estimated parameters were reasonable with the exception of the SC 50 values for ASS and MDH. The scaling factors for all four genes ranged between 1 and 3 and were estimated with good precision.

Class 2
The genes in this group are also up-regulated with acute dosing, but show a biphasic pattern with chronic dosing. The biphasic pattern is more apparent for ornithine decarboxylase (ODC) and CCAAT/enhancer binding protein (CEBP-β) than ASL, suggesting possible existence of dual mechanisms in the regulation of these genes. Model C, where the initial increase in the message is assumed to be mediated by DR(N) followed by a decline driven by a biosignal or a feedback regulator, best described (among other models tested) the general trend of the data for all three genes ( Fig. 6) with relatively good precision for the estimated parameters. Since the biphasic trend is not quite apparent for argininosuccinate lyase (ASL), Models A and B were tested for this gene, however, Model C resulted in best fi t of the data from visual inspection and model selection criteria such as AIC and SC values. The model-fi tted parameters for all three genes are given in Table 3. Since the transduction parameter, k e for CEBP-β could not be estimated with reasonable precision, it was fi xed to an optimum value of 0.1 (upon performing sensitivity analysis).
The degradation rate constant, k dgr,Enzm for ODC (0.14 hr −1 ) was estimated to be somewhat lower than previously reported (0.3 hr −1 ) by Jin et al. (Jin et al. 2003), however, use of a different as well as more complex model could contribute to this difference. Although the IC 50 values for inhibition of the message by the biosignal were quite similar between these genes (ranging from 65-76.6 nM), the values for the amplifi cation factor, γ, were much higher than 1.0 and had signifi cant variation ranging from 3.9 for ASL to 14.7 for ODC. The scaling factors for this group of genes were close to 1.0 with CV less than 20%.

Class 3
This group of fi ve genes showed an early downregulation following acute MPL and a later and sustained up-regulation. However, during chronic infusion stimulation of the messages was predominant, except for glutamate dehydrogenase (GDH) and ornithine transcarbamoylase (OTC).
It is important to note that for chronic infusion, the fi rst sampling was at 6 hr following initiation of infusion. Therefore, it is possible to miss the initial down-regulated phase seen in the acute response profi les. Models D and E, where the nuclear drug-receptor complex (DR(N)) mediates the initial decline followed by an up-regulation caused by a mediator, BS, derived from DR(N), could adequately capture the biphasic patterns of these genes both after acute and chronic administration. A linear stimulation constant (Model D) was suffi cient to describe the biosignal mediated increase in OTC, alanine aminotransferase(ALT) and carbamoyl phosphate synthetase (CPS). However, similar to genes in Class 1 (ASS and MDH), we observed a lack of tolerance in arginase (ARG) and GDH. Model E with saturable stimulation factors (S max and SC 50 ) yielded better fi ttings compared to Model D for these two genes when evaluated by various fi tting criteria. The estimated parameters from model fittings are shown in Table 4. For the three genes which were characterized by Model D, most of the estimated parameters, especially k dgr,Enzm (0.034-1.0 hr −1 ), S DR(N) (0.004-0.027 nM −1 ) and k e (0.03-1.1 hr −1 ) varied between the genes, however, the parameters estimated were usually associated with higher CV%. For the two genes, arginase and GDH (Fig. 8), which could be best described by Model E, the SC 50 parameters could not be estimated with   Table 2. reasonable precision and were fi xed to optimum values based on the various fi tting criteria. Due to the variability associated with data for GDH, both k e and SC 50 values could not be estimated with good precision and therefore were fi xed based on the estimates for the rest of the parameters. Thus the apparent high values for k e (10 hr −1 ) and SC 50 (100 nM −1 ) for GDH should be interpreted with caution. One of the primary reasons for the high CV% of the estimated parameters could be attributed to possible overparameterization of the model with the limited amount of data.

Classes 4A and 4B
Acute dosing of MPL causes down-regulation for this class of genes which can be captured by Model F where DR(N) inhibits the production rate (k syn, Enzm ) of the mRNA. For glucocorticoid receptor (GR) we could observe ( Fig. 9) similar downregulation upon chronic infusion thus acute and chronic dataset could be fi tted simultaneously (Class 4A) using Model F with good fi tting results evaluated by various criteria. For glutaminase (GA) and carbonic anhydrase (CA), although an apparent down-regulation was observed after acute administration without any indication of a biphasic Scaling factor 2.7 (33.9) 1.6 (18.3) 1.14 (34.7) 1.6 (29.6) Enz m0,acute Baseline-Acute 0.89 ( Table 3. pattern, chronic infusion revealed a delayed and sustained up-regulation (Fig. 10) in the gene profi les. This suggests an existence of dual mechanisms after continuous exposure of the drug. No suitable model could explain these differences simultaneously and therefore, only the acute dataset were fi tted with Model F. The estimated parameters for these three genes are given in Table 5. The estimated degradation rate constant for GR mRNA was very similar to previously reported values based on data from Quantitative Northern Hybridization measurements (Sun et al. 1998a;Sun et al. 1998b). Most of the parameters could be estimated with reasonable precision with CV% less than 30%.

Discussion
Corticosteroids, either alone or in conjunction with transcription factors/hormones affect a wide array of genes responsible for regulating various physiological processes. Recently oligonucleotide  Table 4. microarrays were successfully used in our lab to examine these diverse effects of corticosteroids from two sets of animal studies, which allowed us to evaluate changes in an enormous number of genes in various tissues simultaneously after acute (Almon et al. 2004;Almon et al. 2005c;Jin et al. 2003) and chronic MPL. While such evaluations are quite helpful in determining which genes are affected by steroids in rats, another important aspect of such experiments is specifi c disease/   Table 4.
physiological process based evaluation of the data as demonstrated by Almon et al. who examined the effects of MPL on insulin resistance specifi c genes (Almon et al. 2005b) in rat muscle. In this report, we investigated the dynamics of various genes related to the greater urea cycle, an important biochemical pathway present in mammalian livers for nitrogen disposal. Many of these enzymes are known to be highly regulated by CS in the liver (Christowitz et al. 1981;de Groot et al. 1984;Gautier et al. 1985).
The dynamics of many of these genes following acute dosing were previously reported from our lab (Jin et al. 2003). However simultaneous analyses of both the acute and chronic dataset revealed more complex gene regulations which were not apparent after only single high-dose MPL.
Class 1 and Class 4A were the simplest regulation patterns apparent amongst all the classes. TAT, which facilitates the transfer of an amino group from amino acid breakdown to the urea cycle and is also one of the most highly studied biomarkers, could be classifi ed as one of the Class 1 genes. Both TAT and glucocorticoid receptor (Class 4A) dynamics had extensively been studied in our lab (both at mRNA and protein levels) and fi ve generations of highly mechanistic, quantitative PK/PD models were developed for MPL at various dose levels (Ramakrishnan et al. 2002;Sun et al. 1998a;Figure 9. Fitting results of gene array data for glucocorticoid receptor. Solid circles are the mean data with standard deviations. Solid lines are fi ttings with Model F (Eq. (16)) after MPL administration. The estimated parameters are listed in Table 5.  Table 5. Sun et al. 1998b). Thus the presence of these two genes in our database not only had importance with respect to their involvement in the urea cycle, but also allowed us to compare the results from the microarrays to previous measurements by quantitative Northern Hybridization. For TAT the correlations (R) between these two techniques were quite high for both acute (0.68) and chronic (0.87) dosing while for GR the correlation was less than 0.5 which could be attributed to the much lower abundance of this message, the presence of greater variability in these data, and the more limited downward range of observed changes.
The other Class 1 genes were aspAT, with a similar function as TAT, ASS, one of the fi ve urea cycle enzymes and MDH, a citric acid cycle enzyme that catalyzes the conversion of malate into oxaloacetate, thus facilitating recycling of the urea cycle by-product fumerate, which are shown in the urea cycle general schematic in Figure 1. The dynamic profi le for aspAT (Fig. 5), known to be regulated by CS (Garlatti et al. 1996;Pave-Preux et al. 1988) resembled TAT, however, incorporation of a capacity-limited stimulation by DR(N) seemed to provide the best fi t for the data.
The action of CS on ASS and MDH are known to be secondary (Husson et al. 2003;Recupero et al. 1986;Ulbright and Snodgrass, 1993). Bourgeois et al. reported dose-dependent increase in ASS mRNA expression in cultured rat hepatocytes which could be completely abolished by actinomycin D indicating involvement of the synthesis of some secondary factors in this regulation (Bourgeois et al. 1997). The MDH has also been reported to be secondarily regulated by CS by potentiating the binding capacity of T3 receptors since CS action on MDH in thyroidectomized animals was minimal (Recupero et al. 1983;Recupero et al. 1986). However, our current in vivo datasets precluded addition of such secondary factors and therefore for simplicity Model B was used to capture both these enzyme inductions. For both ASS and MDH, however, the general patterns of the stimulation were quite different than TAT or aspAT. While the acute time profi les for ASS and MDH yielded a somewhat delayed and prolonged stimulation (lower degradation rates, 0.15 and 0.08 hr −1 compared to TAT and aspAT, 0.506 and 0.81 hr −1 ), chronic infusion profi les for these two genes did not follow the apparent tolerance phenomenon visible in the chronic DR(N) profi le (Fig. 4). Rather the message seemed to remain up-regulated during the entire period of seven days (Fig. 5) indicating somewhat higher sensitivity to DR(N) compared to aspAT or TAT. The estimated values of the SC 50 parameters for these two genes (Model B) agreed with this observation, almost 15-20 fold (3.2-3.5 nM) lower than the steady-state DR(N) after chronic MPL (57 nM) compared to the higher estimated value for aspAT (191.5 nM).
Reports in the literature show CS-mediated direct or indirect increases in all three genes of Class 2, i.e. ASL, CEBP-β and ODC Hirvonen et al. 1988;Matsuno et al. 1996;Nebes and Morris, 1988). Although our acute dataset predominantly showed up-regulation for these genes, chronic profi les indicated biphasic responses, especially for ODC and CEBP-beta. The initial increase in these messages was followed by a down-regulation up to or close to their baseline followed by a delayed increase. This was similar to the Class 1 genes suggesting a possible existence of a negative feedback regulator limiting the degree of increase of the message during chronic infusion of MPL. Models with a secondary mediator arising from the message itself regulating the initial increase would probably be most mechanistic; however, model convergence could not be achieved with such a model. ASL, one of the fi ve urea-cycle enzymes known to be secondarily regulated by CS, showed similar profi le as TAT or aspAT, however incorporation of the biosignal (Model C) was necessary for satisfactory fi tting of the data. An important transcription factor CEBP-β, contains a GRE in its promoter region and has been postulated to be one of the primary mediators of glucocorticoid effects on urea cycle enzymes in the liver, especially arginase and CPS. This chain of effects, i.e. of stimulation of CEBP-β by CS and in turn delayed stimulation of ARG and CPS by CEBP-β could be nicely captured for the acute dataset (results not shown); however, this model could not be implemented for the chronic dataset where the regulation seemed to be more complex, perhaps due to the involvement of multiple regulatory factors.
Class 3 contains the largest number of genes among all the classes (including three of the fi ve main urea cycle enzymes) predominantly showing a biphasic response, usually a rapid downregulation followed by an up-regulation, which could generally be described by Model D or E. Such a phenomenon is not surprising for CS since they may induce or inhibit the same gene depending on other regulatory proteins/transcription factors that are already present on the DNA control regions (Alberts et al. 2001). Since the upregulation was somewhat delayed for all the genes, both models assume the increase is mediated by a biosignal.
Although most of the reports in the literature suggest an inhibitory effect of CS on OTC in fetal as well as in adult rats (Gautier et al. 1985;Gautier et al. 1977;Ulbright and Snodgrass, 1993), we observed a subtle and delayed up-regulation in this gene after both acute and chronic MPL. As has been mentioned in the previous section, the coordinated regulation of CS on CPS (the fi rst enzyme of the ornithine cycle, converts ammonia, originating from amino acid degradation, into urea) or arginase (the last enzyme of the cycle) via CEBP-β could not be applied to the chronic dataset and thus these two genes were modeled independently of CEBP-β. Arginase and GDH seemed to show saturability in stimulation by the biosignal similar to a few of Class 1 genes. Although this phenomenon could be described by Model E, the precision of the estimated parameters was poor perhaps because of over-parameterization of the model and inadequate data thus limiting the capability of the estimation, especially the capacity-limited parameters such as SC 50 .
Class 4B genes, unlike Class 3 genes, do not show any delayed stimulation after acute MPL, however chronic dosing predominantly showed delayed up-regulation (glutaminase, involved in interconversion of two amino acids, glutamine and glutamate) or a biphasic response (carbonic anhydrase, responsible to form the bicarbonate ion from carbon dioxide from the Kreb's cycle for the fi rst step of urea cycle) suggesting differential regulation between acute and chronic dosing. Not much information is available regarding CS regulation of these two genes in the literature; however, initial investigation of the acute dataset showing downregulation only was surprising to us owing to the fact that CS treatment seemed to increase most of the genes related to the urea cycle either instantaneously or with a delay, which could be clearly observed only with chronic dosing for these two genes.
An important component of our analysis of both acute and chronic dosing of MPL was the necessity of including a scaling factor for the chronic dataset. Since this factor was only included in the chronic dataset (refer to Eq. 17), a value of greater than 1.0 signifi ed greater sensitivity for the chips for the chronic dataset. Two different chips were used for the two treatment groups and data generated from them may not always produce concordant results. In general signal values tend to be higher for the majority of RAE-230A probe sets relative to the corresponding probe sets on RG-U34 suggesting a need for a scaling factor to accommodate this discordance. It was suggested by Affymetrix that the newer 230A chip is most likely to outperform the RG-U34 probe set when both magnitude of signal and responsiveness to a biologically diverse tissue panel are evaluated. The 230A chips were superior in 51% of the probe sets, 5.7% of the probes seemed to have higher signals in RG-U34 chips, and less than 50% of the probes had equivalent signals.
Consistent with the physiological role of glucocorticoids in stimulating gluconeogenesis/urea cycle activity, most urea cycle genes exhibited a general enhanced expression following CS treatment. However, this rich time series data set illustrates that examination of changes in gene expression at a single time following drug administration can be misleading as to magnitude or even direction of change. For example, if a gene such as OTC which exhibits a complex pattern is examined at an early time following acute administration, one would conclude down-regulation by CS. On the other hand, if later time points were examined upregulation is evident. This study also illustrates that time of exposure to drug (i.e. acute versus chronic treatment) is an important consideration to interpretation of the effects of CS on gene expression. Chronic exposure to CS in some cases results in more complex biphasic patterns of expression not evident following acute exposure, indicating probable secondary factors in control of expression. PD modeling can reveal insights into the complexities of mechanisms of gene regulation beyond simple up-or down-regulation by the drug.
In conclusion, endogenous and exogenous corticosteroids affect many physiological processes. Studying each pathway individually is more rigorous and quantitative using techniques such as Quantitative Real-Time PCR, however, it is more time consuming and laborious. Although challenges in proper analysis and understanding of microarray results have emerged since the technology evolved rapidly over the years, we were able to use this plethora of data to understand and analyze the underlying mechanism of various biological processes by using our mechanism-based PD models. The proposed models, while often premised on obvious mechanisms and supported by literature fi ndings, also present testable hypotheses that can be examined in more extensive studies of specifi c genes. Note * Supported by Grant GM24211 from the National Institute of Health